UW-ThPh-2012-17 
AEI-2012-36 



Evolution of a periodic eight-black-hole lattice in 
numerical relativity 

Eloisa Bentivegna 

Max-Planck-Institut fiir Gravitationsphysik 

Albert-Einstein-Institut 

Am Muhlenberg 1, D-14476 Golm 

Germany 



Mikolaj Korzynski 

Gravitational Physics 
Faculty of Physics 

University of Vienna, A-1090 Vienna 
Austria 

Abstract. The idea of black-hole lattices as models for the large-scale structure 
of the universe has been under scrutiny for several decades, and some of the 
properties of these systems have been elucidated recently in the context of the 
problem of cosmological backreaction. The complete, three-dimensional and 
fully relativistic evolution of these system has, however, never been tackled. 
We explicitly construct the first of these solutions by numerically integrating 
Einstein's equation in the case of an eight-black-hole lattice with the topology of 
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1. Introduction 

The vast amount of large-scale cosmological data collected in recent decades has 
shaped a generally coherent picture of our universe [I], where the thermodynamics 
and nucleosynthesis in the early universe, the generation of the seeds of cosmic 
structure and their subsequent evolution all fit together in a simple framework based 
on remarkably few principles. This success, however, comes at the price of accepting 
the existence of a dark sector, i.e. two fluids, dark matter [2j and dark energy |3l|4], 
which have rather peculiar physical properties, have an uncertain collocation within 
the Standard Model, and have never been observed in terrestrial laboratories despite 
accounting for over 95% of the energy density of the universe. 

Whilst this result could very well point to the existence of new physical 
constituents or principles, the possibility that the current way we model the 
inhomogeneous universe be too elementary (a possibility that was advanced for the 
first time in [SI) has now resurfaced [B]. In particular, the question of the extent to 
which cosmic inhomogeneities can dress the value of the cosmological parameters is 
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under scrutiny in a variety of approaches [II [3 El HOl HH [El [13] (see also the review 
articles [UllTn] and references therein). 

An interesting class of models that has been studied for some time is that of 
regular lattices of black holes [IS1[IZ1||] These representations of our universe avoid the 
issues related to the behavior of relativistic fluids (and, in particular, the corresponding 
singularities) ; one could argue that they also represent a more realistic picture of our 
universe, composed of a collection of pointlike objects surrounded by vacuum rather 
than a homogeneous and isotropic fluid with small-scale perturbations. 

Collections of black holes also have the added benefit to be one of numerical 
relativity's routine application areas |18j . from which formalisms, tools and practices 
can be readily imported. In this work, we construct the initial data and simulate the 
evolution of a special sort of black hole lattices, those with extrinsic curvature that 
is initially zero. In section [2] we illustrate how to construct exact initial data for a 
generic black-hole lattice based on the usual Lichnerowicz-York framework (THl [2D] . 
We then discuss a coordinate transformation that simplifies the numerical treatment 



in section 2^ illustrate the details of the evolution in section l3l interpret the results 
and contrast them to the homogeneous and isotropic class in[4j and finally conclude 
in section [Sj Unless otherwise stated, greek indices run from to 4, latin indices run 
from 1 to 3, and we set G — c ~ 1. 



2. Constructing a periodic black-hole lattice 

Given the standard 3-1-1 split of the metric tensor into the spatial metric 7^- and 
extrinsic curvature Kij, initial data for the gravitational field can be generated by 
solving the Hamiltonian and momentum constraints, which read: 

R + K^ - K,jK'^ = IGirGp (1) 

DjKl - DjK = S^Gji (2) 

where R is the scalar curvature of the spatial metric and Di represents the covariant 
derivative associated with 7.^ ; p — n'^n'^T^jy and /' = n^Tj^j^ represent the energy 
and momentum density respectively. 

Let us assume that ji vanishes. A powerful scheme to generate solutions of this 
system is the conformal transverse-traceless framework, which entails the introduction 
of a conformal transformation in the spatial metric, along with the separation of the 
extrinsic curvature into its trace K and traceless part Aij-. 

Itj = V'"* It] (3) 
K 

= y 7»j + (4) 
In terms of these, the constraints take the form: 



R . 1 7 

12 



- ^ ^ - + ^A^.A^^i^-' = -2^G V' p (5) 



b.A'i - ^V'f ■'A/^ - (6) 

X Notice that the usual definition of a black hole in an asymptotically-flat spacetime is inapplicable 
to these spaces. Here, by black hole we denote a spacetime region surrounded by a marginally 
outer-trapped tube (MOTT). 
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Figure 1. The elementary cell of a periodic 3-space containing a number of 
punctures and corresponding inner boundaries, e.g. in the shape of spherical 
surfaces Si. 

A being the laplacian operator of the conformal metric jij, and Aij being related to 
A,j by A,j = ij'^A.j g 

Let us focus on the Hamiltonian constraint. We would like to solve this equation 
with periodic boundary conditions. We also allow for matter content in form of 
ordinary matter p as well as in "punctures", i.e. singularities in ijj of the form mi/2r, 
rui > 0. It can be easily proven that, unlike in the asymptotically-flat case, if we set 
both Kij and R to zero, then this is a slice of Minkowski spacetime. 

To see this, let us first integrate both sides of equation ^ over the fundamental 
cell D (which, for the sake of illustration, we will assume cubical) of the desired lattice, 
with small balls around the punctures excised at the surfaces Si (see Figure [T]). The 
volume integral of A^/' can then be turned into a surface integral on the Si alone, as 
Kij) identically vanishes on the periodic boundary. We thus obtain: 

^ (^f ^ + ^ - A^' A,)j V^d^x = 

§ A different scaling of p is preferable if we want to solve the initial value problem for some types of 
matter, but this plays no role in our argument. 
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= 27rG y^d^ai + S^^im^j . (7) 

On the right hand side, which is manifestly positive, we have the total energy content 
of the cell both in form of a continuous distribution as well as in the punctures. If 
we set both R and K to zero, the equation becomes impossible to satisfy unless the 
matter content vanishes as well. Thus, for non-zero m^, we need to admit either a 
non-zero extrinsic curvature K or a positive spatial scalar curvature R. 

In this work we concentrate on the second case, setting K to vanish. In this 
case the momentum constraint is trivally satisfied and the Hamiltonian constraint 
remains linear in tp, which allows for constructing multiple-black- hole solutions by 
superposition. 

We would like the conformal metric 7 to be periodic just like the physical one. 
The simplest way to ensure that it admits a discrete symmetry is to assume that it is 
a hyperbolic, spherical or flat metric. Since the integral of R must be positive, 7 must 
be the metric of a round 3-sphere. This condition limits both the form of the metric 
tensor and the topology of the lattice [jj] 

Let us note that a similar reasoning applies to the Friedmann-Lemaitre- 
Robertson- Walker (FLRW) models: the Hamiltonian constraint in this class reads 

R 

8 + 12 - 2-^^- 

Since the RHS, representing the matter content, is manifestly positive, then either we 
must have non- vanishing Hubble parameter K, or a positive curvature R (or both). 
Thus, if we assume K — 0, then there is a similarly strong restriction on the metric 
tensor and the topology of the constant-time slices. 

2.1. Punctures on 

Following [5T] we consider puncture-like solutions of the Hamiltonian constraint 
when 'jij and R are respectively the metric tensor 7^^- and the scalar curvature of the 
round 3-sphere: 

(9) 

We fix coordinates on S'^ such that: 

7^^ = dA^ -I- sin^ A {d9^ + sin^ dip^) (10) 

Let us imagine that the sphere is embedded in with the equation {X^)^ -f {X'^)^ -\- 

{X^^^ -\- {X"^^^ = 1; a bar over a capital letter denotes a vector in this space. 

Equation ([9| has no regular solutions, but it is straightforward to find its solutions 
with a puncture-type singularity: 

V'(A)-^479- (11) 
sm A/2 

II In a more general setting, the scalar curvature does not have a definite sign and the integral 
condition ([7| gives little information. Nevertheless, thanks to the Yamabe theorem, 7 must be 
conformally equivalent to a constant-curvature metric. If we rewrite equation ([7| in terms of this 
metric, it obviously becomes a condition for the sign of the Yamabe energy £(7). 
^ The analysis of this initial-data construction in the context of the backreaction problem has also 
recently appeared in I22| 
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or in Cartesian coordinates 



iV=(0,0,0,l) (13) 
We can easily superimpose N such punctures centered at chosen locations 



The parameters > measure the singular part of the solution at the points Ni: 
the leading part behaves Ukc 2Ai/\i. 

Notice that, if one seeks only the regular arrangements of black holes on S^, there 
are only six possible values of A'', corresponding to the six regular tessellations of the 
3-sphere: iV = 5, 8, 16, 24, 120, 600. In the following, we concentrate on the 8-vertex, 
16-cell solution, where the pimcture locations are given by; 

iVi = (1,0,0,0), 
N2= (-1,0,0,0), 
TVs = (0,1,0,0), 

7V4 = (0,-1,0,0), 
iV5= (0,0,1,0), 
7V6= (0,0,-1,0), 
iV7 = (0,0,0,1), 

TVs = (0,0,0,-1). (15) 

and all Ai = 1. The configuration obviously has the symmetry of a 16-cell. In 

particular, it has a discrete group of symmetries generated by tt/2 rotations around 
all pairs of axes of R** and reflections about all four hyperplanes perpendicular to the 
axes. The elementary cell in this pattern is cubical in shape, i.e. it has 6 faces, 8 edges 
and 8 vertices, at which exactly 4 edges meet. All edges lie on great circles of and 
their length is equal to 168.343. 

2.2. Stereographic projection of 

Since it is easier to perform the evolution of asymptotically-flat data as opposed to 
data on a sphere, we employ the stereographic projection from the top of the sphere 
into R^, given by 

2 

= (16) 

It is well known that the projection is a conformal mapping in the sense that 

7f^ = (|f|V4+l)-'5,^ (17) 



or 



2/ Vl-cosA 
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where Sij is the flat metric. The physical metric (14) projected down to R"^ takes the 
form of 



4>{x) 



^ 2A,^l + \n,\y4 



i=2 



\X - Ui 



(19) 
(20) 



Thus the potential consist of iV — 1 punctures of 1/r type, one of the punctures having 
been projected out to infinity. Note that the physical metric involves the (scale-setting) 
factor (Ai)^. We can absorb it by introducing new, rescaled coordinates y — {AiY x. 
The projected conformal factor takes now the form of 



1 



i=2 

N 



=2 2 y-M 



(21) 
(22) 



with rescaled positions of the punctures A^i = {AiY fii. The mass parameters of the 
punctures take the form of 



(23) 



They have the dimension of mass, but do not correspond exactly to the ADM mass 
of the individual punctures measured at their infinities. For the black hole at Ni , for 
instance, the ADM mass is equal to 



N 



N 



^ADM 



IMI 



^ = 4Ai ^ W 1 

i=2 i=2 ' 

This can also be expressed in terms of the original solution tp: 



N 



^ADM ^ 



l-N,-Ni 



(24) 



(25) 



Analogous relations hold for other black holes. 

If we project the 8- vertex solution ( 15 1 down to R'^, it becomes an asymptotically- 
fiat configuration described by ( 22 ) with 7 punctures at points 



A^2 = (0, 0, 0) , 


(26) 


A/3 = (2,0,0), 


(27) 


AA4 = (-2,0,0), 


(28) 


- (0, 2, 0) , 


(29) 


A/g = (0,-2,0), 


(30) 


A^7 = (0,0,2), 


(31) 


A/g = (0,0,-2). 


(32) 



and with mass parameters m2 = 4 and m^, . . . , my — 4\/2. The vertices, edges and 
marginally outer-trapped surfaces (MOTSs) projected to R^^ are presented on Figure[2] 
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Figure 2. Initial MOTSs and elementary cells for the 8-black-hole configuration, 
projected to R'^. The marginal surface corresponding to the black hole at infinity 
encompasses the whole configuration. Note that the 8 cubical lattice cells are 
isometric after the conformal rescaling. 
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3. Evolution 

We perform the three-dimensional evolution of this initial-data set using the Einstein 
Toolkit [23]; in particular, we use the McLachlan code US] to perform a 
finite-difference evolution of the Einstein's equation with adaptive-mesh-rcfinement 
capabilities provided by Carpet [26j. We have also made use of AHFinderDirect [37] 
to search for trapped surfaces. 

Details of the evolution scheme can be found in [Appendix A[ the numerical error 
analysis is performed in |Appendix B[ We evolve the initial data presented in section [2?2| 
in a cubic box of side 40M (here and in the following, we will adopt M = m2/4 as 
the unit of mass, length and time), with a spacing Aq = IM. We refine the grid at 
the seven black-hole locations using eight concentric grids for each location, down to 
a resolution of Ag = = 0.00390625M. 

The eight initial MOTSs have been shown in Figure 2\'^ Notice that, due to the 



stereographic projection, the large marginal surface initially at a radius of 20M is inner 
trapped, rather than outer trapped, i.e. it is the expansion of the ingoing null normal 
that vanishes while the expansion of the outer normal is positive. It is therefore not 
a common apparent horizon of the sort usually encountered in binary mergers. This 
has an interesting side effect: the outer boundary conditions are causally disconnected 
from the region enclosed by this surface. 

We find that the MOTSs at (±2, 0, 0), (0, ±2, 0) and (0, 0, ±2) take approximately 
130M (in coordinate time) to merge to the MOTS at the origin, and approximately 
170M to merge to the larger, inner-trapped one, as illustrated in Figure |4] (the 
asymmetry here is due to the non- uniform numerical slicing). The evolution of 
the z — sections of the marginal surfaces are shown in Figure |3| while the mean 
coordinate radii and masses for the the black holes initially at (0,0,0), (2,2,2) and 
infinity are plotted in Figures [5] and |6] 

From the cosmological standpoint, we are particularly interested in the scaling of 
lengths with proper time. There are at least two candidate quantities for measuring 
the scaling of distances in this system: the (minimal) proper distance between near- 
neighbour surfaces and the proper length of each cell's edges. Note that these estimates 
need not agree with each other, since the expansion rate may well be different at 
different points. 

In order to calculate proper lengths as functions of proper time, we restrict our 
attention to the 1-1-1 subspace spanned by a reprentative curve in time, and obtain the 
proper time and the x-coordinate of gaussian observers with the following relations: 

T{t,x) = [ a{t\x)dt' (33) 



JO 

Xg{t, Xinit) Xg{t ~ At, Xinit) " / P"" {f , Xg{t - At, a;i„it))di' (34) 

Jt-At 

where Xinit is the location of the observer at t = 0. 

The proper distance between marginal surfaces, as a function of proper time, is 
then given by: 

d{t) = / [{-a^T,e) + p^{T,e)){ditf + i3,{Tj)ditdix' + j^,iT,£)dix'dixq^^^di 

J It 

+ We really only track the surfaces corresponding to the black holes at the origin, at infinity, and on 
the positive x-axis. The locations and shapes of the remaining five surfaces, included in the plots for 
clarity, are obtained by symmetry arguments. 
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Figure 3. Section of the eight marginal surfaces on the 2 = plane, at times 
t = (0, 25, 50, 75, 100, 125, 150)M. The bottom plot is a zoomed-in version of the 
seven central MOTSs only. 



Evolution of a periodic eight-black-hole lattice in numerical relativity 



10 



0.25 




-0.25 I ' ' ' ' ' ' ' ' ' 1 

-0.25-0.2-0.15-0.1-0.05 0.05 0.1 0.15 0.2 0.25 

x/M 




Figure 4. Section of the eight marginal surfaces on the z = plane, on the spatial 
slice of the first surface merger, corresponding to coordinate time t = 128M (top) 
and on that of the second merger at t = 167M (bottom). 
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Figure 5. Masses of the three marginal surfaxies of the black holes at (0,0,0), 
(2,2,2) and infinity, as a function of coordinate time. 
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Figure 6. Mean coordinate radii of the three marginal surfaces of the black holes 
at (0,0,0), (2,2,2) and infinity, as a function of coordinate time. 
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Figure 7. Several measures of scaling in the eight-black-hole universe, as 
functions of proper time t, plotted against a possible identification of the 
corresponding FLRW model (see section |4] for details). All the quantities have 
been renormalized to their respective values at r = 0. 



(35) 

where 77- is the shortest constant-r geodesic, parametrized by £, connecting two 
surfaces. We measure this quantity for the two outer-trapped surfaces initially at 
(0,0,0) and (2,0,0), the geodesic lying on the x-asds for symmetry reasons. This 
quantity is plotted in Figure [7] 

Similarly, the proper length of a lattice edge is given by equation ( 35 ) , where now 
7t- is the constant-r geodesic, parametrized by £, connecting the two vertices. It is easy 
to work out the initial shape of the cells, illustrated in Figure [2j The initial locations 
of the 16 vertices are given by (±2/3, ±2/3, ±2/3) and (±2, ±2, ±2). For simplicity, we 
choose to focus on the edge connecting (2/3,2/3,2/3) to (2,2,2), which, for symmetry 
reasons, always lies along the x — y — z diagonal. The edge's proper length as a 
function of proper time is also shown in Figure [Tj 

For reference, the relative spatial and temporal scales of the system are illustrated 
in Figures [8] and |9] The numerical locations of the two vertices during the evolution 
is shown in Figure [8j along with a few other representative points on the geodesic and 
the constant-r lines. In Figure [9j we also show the span of the numerical coordinates 
of the cell edge in (r, Xg) space: this illustrates how the gauge condition adopted in 
this simulation freezes the evolution around r < 150M, preventing us from observing 
the system's behavior after this time. In Figure |8j we also plot the intersection of 
the marginal surfaces surrounding the black holes at the origin and at infinity with 
the X = y = z diagonal; this illustrates that, after t ~ 120M, the inner vertex is 
quite close to the MOTS of the black hole at the origin: by this time, we can expect 
finite-size effects to play a significant role in the scaling of lengths. Additionally, we 
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0.0001 



Figure 8. Coordinate Xg of gaussian observers and constant-r lines (r = 
0, 10, . . . , 140M) on the x = y = z diagonal, for points initially located between 
the two vertices at (2/3,2/3,2/3) and (2,2,2). The thick lines represent the x- 
coordinate of the intersections of the marginal surfaces at (0,0,0) (left) and at 
infinity (right) with the diagonal, up until the mergers at t = 128M and t = 167Af, 
respectively. 



show in Appendix B tliat the numerical error quickly degrades after r ^ 80M. Based 
on these considerations, we only show the scaling up to proper times of about lOOM 
in Figure [Tj 

Figure [7] also includes a possible counterparts of the eight-black-hole lattice in 
the FLRW class: one with the same initial edge length. This definition is made more 
precise in section [4j where we will constrast these results with the large-scale dynamics 
of a perfectly homogeneous and isotropic universe. 

Finally, let us notice here in passing that the interest of this toy model goes 
beyond the cosmological application. In particular, it provides an interesting example 
of overlapping MOTSs within the framework of a BSSN evolution. 



4. Comparison with the FLRW class 

The comparison of the configuration with an FLRW model requires solving Ellis' 
"fitting problem" i.e. determining the parameters characterizing the reference 
FLRW model which our configuration resembles most closely. There are infinitely 
many ways to do this; we will sketch below the procedure we will use in this paper 
based on the quantities measured in section [3] 

Due to the symmetry group of our configuration, the reference FLRW model is 
a closed model (fc = 1), with spatial slices of spherical topology. The matter content 
is represented by dust. Since the primary variable describing any FLRW model is the 
scale factor, representing the scaling of lengths in time, the first step in the comparison 
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Figure 9. Space-time region spanned by tlie edge between (2/3,2/3,2/3) and 
(2,2,2), represented in (r, Xg) coordinates. The slicing adopted in this work only 
extends, in this region, up to r < 150M. The horizontal lines are constant-t lines, 
i.e. sections of the spatial hypersurfaces used in the simulation. 



is to identify some measure of length in the lattice universe. We required this variable 
to be non-local, thereby capturing the large-scale, average behaviour of the universe 
rather than the local physics at a single, arbitrarily chosen point. The total volume 
of the configuration, which is the most obvious parameter for a closed universe, is 
obviously infinite and therefore of no use for our purpose. 

The problem of the size measurement is a bit simplified by the discrete symmetry 
of the model. In section [3] it seemed reasonable to choose the variable in a way 
which is consistent with the cell structure of the black hole lattice. The first obvious 
choice would be to use the geodesic distance -Dhor(T) between the MOTSs of two 
neighbouring black holes. Recall that a MOTS is a closed two-surface whose null 
expansion vanishes in one direction [28 . We have observed that the behavior of this 
quantity for small times varies considerably from the behavior of the size parameter 
of a closed FLRW (see Figure [?]). In particular, the derivative of -DhorlT) does not 
vanish at t = 0, but rather approaches the limiting value -Dhor(O) = —2, despite the 
configuration being momentarily static and symmetric with respect to time reflections. 
This striking phenomenon can be explained by the fact that the black hole MOTSs at 
r = are all bifurcation surfaces of the horizon. One can check that in our coordinates 
the corresponding MOTSs seem to approach each other at approximately the speed 
of light even at the moment of maximal expansion (see Appendix C). This makes 
-Dhor(T) an unsuitable size parameter for the purpose of FLRW fitting. 

Our second choice for the size variable was the geodesic length of the individual 
cell's edge, i.e. the geodesic distance between the two vertices of an individual cube, 
also used in The edges lie relatively far from the punctures, and thus from the 
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potentially problematic black hole region, at all times. Hence we can hope they are 
resolved quite well during the simulation. 

The mapping to the FLRW class is then carried out by fitting the size acs{t) of 
the reference sphere discussed in section [O] by demanding that the length of the edge 
match exactly the length of the corresponding edge of the cubical tiling of a round 
(FLRW) sphere. At the same time this fitting gives the effective Ricci curvature '^^^Rcs 
as the curvature of the corresponding FLRW sphere: 

(3)i?,ff = A- (36) 

At first sight, it may seem quite strange that the effective Ricci curvature is not 
given by any kind of average, over a domain in the constant time slice, of the local 
values of the curvature. Keep in mind however that in the FLRW class, due to the high 
symmetry, the value of the (constant) Ricci curvature is directly related to infinitely 
many other parameters characterizing the geometry, for example the relation between 
the volume and the area of spheres, the angle deficits, the geodesic focusing at any 
point, and so on. Since the geometry of the discrete models is not homogeneous, 
these simple relations are lost. Nevertheless it is not a priori clear which of those 
parameters we should regard as a convenient inhomogeneous generalization of the 
spatial ^^^R appearing in the Friedmann equation. The definition of ^^^Rcs we propose 
here is based on rescaling the value of the Ricci curvature of a unit 3-sphere by the 
ratio of the size of the lattice edges. As we shall see, it provides, together with 
a complementary definition of the effective energy (see below), an excellent fit for 
the dynamics of the inhomogeneous lattice until up the time of t « 80, when we start 
losing resolution (mostly due to an outgoing shift vector that makes the system rapidly 
shrink in coordinate size). Independently from this limit, after t ~ lOOA/ the eight 
MOTSs rapidly swallow up the whole spatial slices, making this model less and less 
appropriate to describe a universe filled with point-like masses. 

In order to derive an effective energy density, let us recall that the configuration 
is at rest at t = 0, so we can assume that Oeff (0) = and use Friedmann's equation to 
obtain poff 

0=y/3off Q — ■ (37) 

It is instructive to compare the effective total mass, obtained as the product of poff with 
the volume of the FLRW sphere, with the total matter content of our configuration: 

Mefi = poff 2TT^alff = 378.78, A/sbh = 8i\/ADM = 303.53 (38) 

Clearly the effective mass is around 25% larger than the sum of ADM masses of 
the individual black holes. This is consistent with the expected nonlinear effects of 
gravitation, as the strong gravitational fields gravitate themselves. It also agrees quite 
well with (21] , where the fit with the corresponding FLRW model is based on imposing 
the equality of the total masses and the lengths of the edges are used for comparison 
instead. 



5. Conclusions 

We have discussed the construction of the initial data corresponding to periodic lattices 
of black holes using the Lichnerowicz-York construction. We have shown that just like 
in FLRW class there is a link between, on one hand, the curvature of the underlying 
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constant-curvature metric (and thus also the topology of the spatial slice) and, on the 
other hand, the matter content and the Hubble parameter. In particular, it follows that 
it is not possible to find a flat, rectangular lattice of black holes without a momentary 
expansion or contraction. 

We then recalled the construction of periodic lattices of black holes on an 
sphere with vanishing extrinsic curvature, originally introduced in |21j . We focused 
on the symmetric eight-black-hole configuration and showed that it can be conformally 
projected to an asymptotically-flat seven-black-hole configuration with unequal mass 
parameters. We have discussed briefly its basic properties and then showed the results 
of the numerical evolution of this system. The main goal was the comparison with 
the time evolution of a dust-filled closed FLRW model, evolved from the maximal 
expansion moment until the recollapse. 

We proposed two methods for measuring the effective size of the configuration: 
the geodesic distance between the MOTSs of two neighbouring black holes and the 
geodesic length of the edge of a single cell. The first one turned out to behave very 
differently from the size parameter of FLRW, as its derivative did not seem to vanish 
even at the initial slice. We explained this unusual behavior by noting that, at the 
maximum expansion, each MOTS is a bifurcation surface where two distinct marginal 
tubes intersect. The edge length, on the other hand, calculated in normal coordinates 
and in proper time, seems to follow very closely the evolution of a closed FLRW if we fit 
its size (and consequently its curvature and mass) in an appropriate way. Despite our 
configuration being very far from homogeneity, the effective size obeys the Friedmann 
equation to a remarkable degree up to times of t 80M, which is approximately 30% 
of the recollapse time of the FLRW. After that time our simulation is simply unable 
to resolve the system. The only observable backreaction (or coarse-graining) effect in 
our simulation seems to lie in the effective total mass of the system, which turns out 
to be 25% larger that the sum of masses of the individual black holes. In other words, 
the eight-black-hole lattice does mimic a closed FLRW dust model, but one whose 
total mass is substantially larger than that due to the black holes alone. 

It is of course not clear to what extent these results will hold if we consider other 
types of models (flat or open, with a positive initial expansion) or if we drop the 
assumption of existence of a large group of discrete symmetries. Note that at first 
sight the symmetry assumption may look like an innocent ansatz whose only purpose 
is to simplify the geometry of the problem. Nevertheless it is important to understand 
that it is in fact quite restrictive. In particular it prohibits many types of interactions 
between the matter inhomogeneities, such as two-black-hole mergers or interactions 
via \ow-£ spherical harmonic modes. Configurations which do not have that kind of 
symmetry may potentially exhibit many other effects of backreaction. 
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Appendix A. Evolution system 

In order to solve Einstein's equation, we use the McLachlan code, which implements 
a finite-difference discretization of the BSSN formulation ^51 1501 15T] . 

{dt - l3'di)W = -^aK+ ia,/3^ (A.l) 
{dt - p'di)K = - D,D'a + a(A,, A^ + ^K^) (A.2) 

(dt - (3'di)% = - 2aAj + 27zO-5fc)/3* - (A.3) 

[dt - l3'di)A,, = W^{-D,D,a + aR,, f^ (A.4) 
+ a{KA,,-2A,kA'}) (A.5) 

+ 2Afe(,5,)/3'= - ^-A,,dkP^ (A.6) 
[dt - I3'di)r = ^^'^P'djPu + \f'd,dA - t'd.P' (A.7) 
+ \rdjP^ - 2A'^dja (A.8) 

+ 2a(f A^'*^ - SA^'^dk InW- \f^djK) (A.9) 

where 7ij — W^^^ij is the three-metric, Kij — ^^ij + W^^Aij is the extrinsic 
curvature, and W = det7~-'^/^ and = —9^7*^ are auxiliary variables. The gauge 
variables are evolved according to: 

{dt-P'di)a = ^2aK (A. 10) 

idt-(3'd,)l3^ ^^B^ (A.ll) 
{dt - P'di)B= = - {dt- P'd^i)V^ - vB^ (A.12) 

with T] = 1. 

We use fourth-order finite-differencing. An important difference with standard 
black-hole evolutions is that we choose a = 1 everywhere as the initial condition, 



because a precollapsed lapse a = (see (21)) leads to a vast collapsed region at 
the center of the domain, which unnecessarily slows down the proper-time evolution 
of the black holes. 

Appendix B. Numerical error on the proper distance estimate 

Estimating the numerical error associated with a three-dimensional, complex 
simulation is notoriously difficult. Additionally, our scheme for computing the proper 
distance does not only involve finite differencing and AMR operations, but also a 
number of post-processing steps, in particular reslicing 1-1-1 spaces in terms of gaussian 



observers (which involves using (33 1 and (34) to obtain the gaussian coordinates, and 



then an interpolation in both time and space) and integrating the line element to 
obtain a proper distance. 

In order to quantify the cumulative error of this procedure, we evolve the same 
initial data at two additional resolutions, corresponding to a spacing of Aq = 2M 
and A = AM on the coarsest grid, and compare the edge proper length to obtain an 
order of magnitude and a scaling for the numerical error on this quantity. The result 
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Figure Bl. Convergence study for the edge proper length. Dc(t), DmiT), and 
D{{t) are the edge proper lengths for the coarse, medium and fine resolution 
respectively. The factors Fi represent j-th order scaling of the error with the 
lattice spacing. 



is shown in Figure |B1[ the error scales hke the second power of the lattice spacing 
up until t ^ 40M, degrading afterwards. The order of magnitude is below IM by 
t = 80M (or about 0.6% of the proper distance); this seems like a reasonable number 
to represent the numerics-related error bar on the proper distance. 

Appendix C. Geodesic distance between the MOTSs 

In this appendix we will argue that the minimal surfaces around the punctures at i = 
are all bifurcation surfaces, i.e. they give rise to two MOTTs, one propagating towards 
the asymptotically flat end of the manifold and the other in the other direction. Let 
Eq be the t = hypersurface. Consider the evolution of the initial data near t = 
with gauge conditions a = 1 and /3' = at Eq, as we have chosen on our initial 
hypersurface for the numerical evolution. 

We define two null normals to D given by ? = ^ (T + i?) and fc = ^ (T - i?). 



see Figure CI It is easy to see that the expansion of both null normals and 9 
vanishes. The reason for that is the fact that D is minimal, therefore the expansion 
with respect to R vanishes, and that the initial data is time-symmetric {Kij =0), so 
the expansion in T direction is zero as well. 

The behavior of the MOTS is determined by the continuation equation for a 
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Figure C2. Two neighbouring minimal surfaces at t = as bifurcation surfaces. 
The outward- moving MOTTs (solid line) seem to expand with approximately 
light speed. 



MOTT, which gives the condition under which a variation of a MOTS preserves the 
vanishing of one of the nuU expansions, see |321 128j . Let now n = nil -\- nkk he the 
variation vector field. Since both nuU expansions vanish, we may continue the MOTS 
by imposing either the condition 0' = or ^'^ = 0. In the first case the continuation 
equation reads 

A-2lj^Da - (^^+ Dauj^ - UJA ) "fe - ni |' - 0. (C.l) 

A denotes here the Laplace operator on Z?, loa = Vyifc" la-, c; is the shear of / and 7?. 
is the scalar curvature of D. 

It turns out that both the surface D itself and the physical 3-metric 7 in its 
vicinity are almost exactly spherically symmetric. The reason is that D lies relatively 
close to the pucture p, at radius r w 0.21 in the R'^ variables, small comparing to 
the distance of 2 between the punctures. The round 3-sphere metric 7"^ is of course 
spherically symmetric around the puncture. The conformal factor is dominated 
in this region by the singularity at the origin and the sum of six terms coming from 
other punctures. The latter is not exactly spherically symmetric, but the non-spherical 
contributions from those 6 terms cancel out to a great degree due to their symmetric 
alignment. Therefore the relative variation of V' in angular variables at that place is 
only of the order of 10~^. All deviations from spherical symmetry are thus of very 
small order of magnitude, so in the first approximation we may neglect all terms in 



(C.l I which break the spherical symmetry of 7, i.e. those containing the vector field 
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Figure C3. The proper distance between near-neighbour MOTSs as a function 
of proper time. The straight hne has a slope of -2. 



and the tensor field cr/ The resulting equation now reads 

A-|^nfe-0. (C.2) 

The solution is straightforward: = 0, ni being an arbitrary function. This obviously 
corresponds to the continuation of the MOTT along the null vector field /. The same 
argument shows that we may continue the other MOTT along k. Both MOTTs are 
thus in very good approximation isolated horizons. 

It is now clear why the derivative of geodesic distance between MOTSs does not 
vanish despite the configuration being time reflection symmetric. The MOTTs of two 
neighbouring black holes, when viewed in either Gaussian or numerical coordinates, 
seem to expand almost exactly at the light speed (equal to 1 in simulation units). This 
by no means contradicts the time reflection symmetry of the inital data, as the MOTT 
expanding in the future direction is accompanied by another one which expands in 



the past direction, see Figure C2 Figure C3 showing the distance between MOTSs 



along with the expected initial slope of -2, provides support to the analysis above. 



References 

[1] Cervantes-Cota J L and Smoot G 2011 (Preprint [ll07 . 1789^ 

[2] Sanders R 2010 The Dark Matter Problem: A Historical Perspective (Cambridge University 

Press) ISBN 9780521113014 
[3] Wang Y 2010 Dark Energy (Wiley- VCH) ISBN 9783527409419 

[4] Amendola L and Tsujikawa S 2010 Dark Energy: Theory and Observations (Cambridge 
University Press) ISBN 9780521516006 



* Actually the shear term |(T;p in (C.l I, being quadratic, is of even lower order. 



Evolution of a periodic eight-black-hole lattice in numerical relativity 



21 



[5] Ellis G and Stoeger W 1987 Class. Quantum Grav. 4 1697-1729 

[6] Andersson L and Coley A 2011 Class. Quantum Grav. 28 160301 

[7] Korzyriski M 2010 Class. Quantum Grav. 27 105015 

[8] Bolejko K, Celerier M N and Krasinski A 2011 Class. Quantum Grav. 28 164002 

[9] Buchert T 2011 Class. Quantum Grav. 28 164007 

[10] Clifton T 2011 Class. Quantum Grav. 28 164011 

[11] Brannlund J, Hoogen R v d and Coley A 2010 Int.J.Mod.Phys. D19 1915-1923 

[12] Wiltshire D L 2009 Phys. Rev. D 80 123512 

[13] Wiltshire D L 2007 in "Dark Matter in Astroparticle and Particle Physics: Proceedings of the 
6th International Heidelberg Conference", eds H.V. Klapdor- Kleingrothaus and G.F. Lewis, 
(World Scientific, Singapore, 2008) pp 565-596 

[14] Clarkson C, Ellis G, Larena J and Obinna U 2011 Rep. Prog. Phys. 74 112901 

[15] Wiltshire D L 2011 Class. Quantum Grav. 28:164006, 2011 

[16] Lindquist R W and Wheeler J A 1957 Rev. Mod. Phys. 29 432-443 

[17] Clifton T and Ferreira P G 2009 Phys. Rev. D 80 103503 

[18] Pfeiffer H P 2012 {Preprint 1203.5166) 

[19] Lichnerowicz A 1944 J. Math. Pures Appl. 23 

[20] York James W J 1971 Phys. Rev. Lett. 26 1656-1658 

[21] Wheeler J 1983 Foundations of Physics 13(1) 161-173 ISSN 0015-9018 10.1007/BF01889418 

[22] Clifton T, Rosquist K and Tavakol R 2012 (Preprint 1203.6478) 

[23] Loffler F, Faber J, Bentivegna E, Bode T, Diencr P et al. 2011 {Preprint |1111.3344[ | 

[24] Mclachlan code URL http://www.cct.lsu.edu/~eschnett/McLachlan/ 

[25] Kranc code URL http://kranccode.org/ 

[26] Schnetter E, Hawlcy S H and Hawkc I 2004 Class. Quantum Grav. 21 1465-1488 

[27] Thornburg J 2004 Class. Quantum Grav. 21 743-766 

[28] Andersson L, Mars M and Simon W 2007 {Preprint 0704.2889') 

[29] Nakamura T, Oohara K and Kojima Y 1987 Prog. Theor. Phys. Suppl. 90 1-218 

[30] Shibata M and Nakamura T 1995 Phys. Rev. D 52 5428-5444 

[31] Baumgarte T W and Shapiro S L 1999 Phys. Rev. D 59 024007 

[32] Korzyhski M 2006 Phys. Rev. D 74 104029 



